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Abstract 

In this paper we research the extraction of the angular rate vector from attitude 
information without differentiation, in particular from quaternion measurements. We show that 
instead of using a Kalman filter of some kind, it is possible to obtain good rate estimates, suitable 
for spacecraft attitude control loop damping, using simple feedback loops, thereby eliminating 
the need for recurrent covariance computation performed when a Kalman filter is used. This 
considerably simplifies the computations required for rate estimation in gyro-less spacecraft. 
Some interesting qualities of the Kalman filter gain are explored, proven and utilized. 

We examine two kinds of feedback loops, one with varying gain that is proportional to the 
well known Q matrix, which is computed using the measured quaternion, and the other type of 
feedback loop is one with constant coefficients. The latter type includes two kinds; namely, a 
proportional feedback loop, and a proportional-integral feedback loop. The various schemes are 
examined through simulations and their performance is compared. It is shown that all schemes 
are adequate for extracting the angular velocity at an accuracy suitable for control loop damping. 

INTRODUCTION 

In most spacecraft (SC) there is a need to know the angular velocity. Angular velocity is 
needed for two major tasks; namely, attitude control, and attitude computation. When the attitude 
is given from an autonomous star tracker (AST) - for example 1 , then, of course, the need for 
angular velocity exists only for attitude control. Traditionally, the SC angular rate vector is 
obtained from gyroscopes installed on board. In smaller, lighter and low-cost SC, there is often a 
desire to do away with gyros and use other means to determine the SC angular rate. Even in gyro- 
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equipped satellites there is a need to determine the angular rate by other means when the angular 
rate is out of range of the SC gyros or when a total gyro failure occurs. 

There are several ways to obtain the angular rate in a gyro-less SC. When the attitude is 
known, in whatever parameters it is given, one can differentiate it and use the kinematics 
equation that connects the derivative of the attitude with the satellite angular rate in order to 
compute the latter 2, 3 . However, since the attitude is obtained from sensor measurements, the 
differentiation of the attitude introduces a considerable noise component in the computed angular 
rate vector. 

Another approach is that of using the attitude parameters, or the measured directions 
themselves, as measurements in some kind of a Kalman Filter (KF). In this case the kinematics 
equation that connects the attitude parameters, or the directions, with their derivatives, are 
included in the dynamics equation used by the filter, thereby the need for differentiation is 
eliminated 4, 5 . However, the use of a KF of some kind requires the computation of a covariance 
matrix. Not only is this process cumbersome, sometimes it may also pose an accuracy problem. 
These accuracy considerations led to the use of the more computationally intensive covariance 
measurement-update formula 6 , square root filtering 7 , and other sophisticated approaches. 

In this work we investigate simple algorithms that do not require covariance computation. 
In particular we investigate the use of passive feedback loops for extracting the angular velocity 
from attitude information without differentiation. In our approach, the dominant factor is 
simplicity rather than accuracy, because for control loop damping, crude rate information is often 
sufficient. We start our investigation in the next section with an examination of the Pseudo- 
Linear Kalman (PSELIKA) filter for estimating angular rate when quaternion measurements are 
given. In particular, we consider the filter gain when the filter tracks the SC rates well. Some 
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particular characteristics of the gain are observed and analytically proven. In the section that 
follows we introduce a simple realization of the PSELIKA which is based on these particular gain 
characteristics. This leads to very simple passive feedback systems, which are introduced in the 
following section. All the systems are tested with continuous quaternion inputs. Discrete inputs, 
with and without measurement noise are presented in the next section. The tests prove that the 
systems are suitable for rate estimation for control loop damping. Finally, in the last section we 
present our conclusions. 

ANGULAR RATE DETERMINATION BY ESTIMATION 

In this work, we consider the quaternion of rotation as the attitude measurement. The rate 
of change of the quaternion is described by the well-known equation [see e.g. Ref. 8, p.5 12] 

q=|Qq (l.a) 

where 

0 © z -co y co x 
-co z 0 m x co y 

co y -<D X 0 co z 

-®x -® y -® z 0 

and co x , co y , and co z are the components of the angular velocity vector, co , T designates the matrix 
transpose and [-cox] is the cross product matrix of the vector -to . The cross product, [ax] , of 
the general vector a is given by 

0 -a z a y 

[ax] = a z 0 -a x (l.c) 

_ -a y a x 0 

where a x , a y , and a z are the components of a . Equation (l.a) can also be written as 

q=lQ«> (2. a) 

where 


[-cox] co 
-co T 0 
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Q = 


ql 3 + [ex] 


-e 


(2-b) 


The elements q and e are respectively the scalar part and the vector part of the quaternion; that 
is, q T = [e T q] , and I n is the ^-dimensional identity matrix. 

The angular dynamics of a rigid body SC is given in the following equation [8, p. 523] 

I© + h + (0 x (I© + h) = T (3) 

where I is the SC inertia matrix, h is the angular momentum of the momentum wheels, and T is 
the external torque acting on the SC. Since the inertia matrix, I, is invertible we may write this 
equation as 

© = r 1 [(Ito + h)x]to + r 1 (T - h) (4.a) 

f (to) = r 1 [(I© + h)x]to + r 1 (T - h) (4.b) 


Define f (to) as follows 


then we can augment Eqs. (2.a) and (4), and add to them white noise vectors to reflect modeling 
uncertainty. This yields the state space augmented equation 


IQ© 

m 


(5.a) 


In the case considered where the measured entity is the quaternion itself, the measurement 
model is 


q.-[i. “I 


+ v„ 


(5.b) 


where q m is the measurement, and v q is an appropriate measurement noise. Equations (5) have 

been used to estimate the angular rate with much success either after linearization, using an 
extended Kalman filter (EKF) 9,10 or, after a simple manipulation by a pseudo-linear Kalman filter 
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(PSELIKA) 11 . In fact, when the measurements come at a relatively fast rate, satisfactory rate 
estimate can be achieved when the nonlinear SC dynamics model of Eq. (4.a) is replaced by the 
following simple linear Markov model 12 


d) = -Nco + \ a (6) 

in which N = -nl 3 , n is an inverted time constant, and v M is an appropriate noise vector 11 . In 


this case we obtain the following model 


"q" 

__ 

^ 4 x 4 

iQ 

V 

+ 

V 

d> 


_^ 3 x 4 

-nl 3 _ 

to 


TV 


( 7 ) 


The dynamics model of Eq. (7) and the associated measurement model of Eq. (5.b) were used 
successfully for estimating SC angular rates using PSELIKA 11 . PSELIKA is an ordinary linear 
Kalman filter (KF) where state-variables in the dynamics and the measurement matrices, are 
replaced by their current best estimate. 

As seen in the following continuous general linear KF algorithm, 

i = Fx + K[z - Hx] (8) 

a KF of any kind (or a Luenberger observer 13 ) operates as a feedback system. The difference 
between measurement, z , and its estimate, Hx , is multiplied by the gain, and the product is used 
to correct the state estimate rate of change. In our case of rate estimation, when the continuous 
PSELIKA is applied to the dynamics and measurement models that are given, respectively, in 
Eqs. (7) and (5.b), Eq. (8) becomes 

A 

q 

A 

CO 

which can be written as 

It should be noted that this is just a model and not the true description of the rate vector. 



^4x4 

±Q‘ 

* A 

q 

+ 

1 

1 


_^3x4 

-nl 3 

0) 


L K J 


(q»-q) 


( 9 ) 
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( 10 ) 


A block diagram realization of the last equation is shown in Fig. 1 . Note that Q is a function of 
q. In fact, since the AST yields very accurate quaternions, we can use q m rather than q to 
evaluate Q . The Kalman gain components K q and K m are, of course, computed as a part of the 
usual gain computation of the KF. 



Fig. 1: Block Diagram Representation of the PSELIKA Filter Rate Estimator 


Although in practice the discrete formulation of the filter algorithms are implemented, we 
interchangeably use the continuous formulation in order to investigate certain qualities of the 
filters. This is permissible because the effect of a continuous filter can be considered as that of a 
discrete filter running at extremely small step size between measurements. 

Qualities of K q and 

In estimating the angular velocity using quaternion measurements and the PSELIKA filter, 
the gains K q and K m have special forms that could be used to our computational advantage. To 

explore these qualities, it is convenient to consider the continuous KF algorithm as a discrete 
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algorithm. Omitting, for convenience, time labeling, the formula for computing the discrete 
Kalman gain is 


K = P(-)H t [HP(-)H T + R J 1 (1 l.a) 

where P(-) is the a-priori covariance matrix, His the measurement matrix, and Ris the 
covariance matrix of the measurement noise, all at time t k . From Eq. (5.b) 


H=[I, 0] 


For this particular measurement matrix one obtains 


p i,(-)[ p „(-) +R ] 1 

P5(-)[Pi,(-)+R]"‘ 


where P u (-) and P ]2 (-) are sub-matrices of P(-) as seen below. 


(11 .b) 


(ll.c) 


From Eq. (1 1 .c) we obtain 


P(-) 


P.iH P«(-) 

.PnH P 2 2 (-) 


K q =P i I (-)[ P i l ( _ ) + R ] ' 

K m = P? 2 (-)[? u (-) + R]-' 


(lid) 


(He) 

(Ilf) 


Since P u (-) and R are symmetric, from Eq. (1 l.f), it is clear on the outset that K q is symmetric 
too. In order to examine the form that K 0) takes, and to further examine the form of K q , consider 
the transition matrix from time t k _, to t k . Let A = t k - t k _ 1 then A , the transition matrix, is given 
by A = exp(FA) . For small A the first order approximation of the exponential function suffices. 
It is 


A ~ I 7 + 

^4x4 

H— 

o 

i 

A = 

"I 4 

iQA 


_^3x4 

-! 3 n_ 


_^3x4 

I 3 (l-nA) 


(12.a) 
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Assuming that the process noise is white, that the variance of each quaternion state is r| q and 
each angular velocity state is r| (0 , the a-priori covariance matrix at time t, is computed as follows 


P,(-) = A 0 P 0 (+)Aj + 


V4 0 4,3 
L^3 


*q 

^3,4 Pro^ 


(12-b) 


where P 0 (+) is the initial covariance matrix and P,(-) is the a-priori covariance matrix at time 


t, . Following the usual practice, we choose P 0 (+) as 


P 0 (+) = 


Pq,0^4 ^4,3 

^3,4 Pro, 0^3 


Using Eqs. (12), one obtains 


P>(-) = 


(P q +P q )I 4 +7QQ T ProA 2 iQPfflA(l ~nA) 

, iQ T ProA(l-nA) I 3 [(1 nA) 2 P(0 -4- ] 


Neglecting, where appropriate, terms containing A 2 yields 




QP(o A 


(Pq+Pq) 1 . 2 

. IQ T Pro A I 3 [(l-2nA)p.+T|J 


(12.c) 


(12.d) 


(12.e) 


Let R = £I 4 , then using the last equation, K, , the Kalman gain at time t, , is computed as 
follows 


Define 


K, = P,(-)H t [HP,(-)H t + rJ 1 = 


^(Pq+PqVCPq+Pq+S) 

7Q T ProA/(P q +P q +^) 


then K, can be written as 


a = (P q +P q )/(P q +P q H) 
P= 2 Pro /(P q + *!<,+£) 


K,= 


I 4 a 

Q t Pa. 


(12.f) 

(12-g) 

(12.h) 

(12. i) 
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This gain matrix has a very peculiar form. The question is whether this peculiarity is a 
consequence of the fact that P 0 (+) is a diagonal matrix. To answer this, we go through one more 


iteration. To do that, we compute, P,(+) , the updated covariance matrix as follows 

P, (+) = (I- K,H)P, (-)(I - K,H) T + K, • R • K[ = 
I 4 [(l-a) 2 (p q +q q ) q +^ a2 ] 

Q T (ti P<o A - P A (P q + P q )]0 “ a ) + a £P A ) 

Q ([tP<o A - P A (P q + P q )10 - a) + °^P A ) 

I 3 [(l - 2nA)p (0 + - Pp m A 2 + p 2 A 2 ( Pq + q q ) + ^p 2 A 2 ] 


Define the following constants 

y = (1 -a) 2 (p q +p q ) q + £,a 2 


X = [tP<o A - P A (P q + P q )10 “ «) + a^p A 
8 = 0- 2nA)p ro + r| ffl - Pp m A 2 + p 2 A 2 (p q + r| q ) + ^p 2 A 2 


«(l-2nA)p. + ri. 
then Eq. (12.j) can be written as 


P,(+) = 


hi 

Q T x 


Qx 

I 3 5 


(12.k) 


(12.0 


Next we propagate this covariance matrix to time point t 2 
P 1 (-)P 1 (-) = A,P,(+)A, T +0 


■ 1 . iQA 

'I 4 Y 

Qx~ 

I 4 o 4 , 3 

_L 

^4,3 

_0„ (l-nA)I, 

Q T x 

I 3 5 

_iQ T A (l-nA)I 3 

T 

_ ^3,4 


This results in 


P 2 (-) = 


I 4 Y + QQ Ta X + 1 QQ T a2 8 + rjql 4 
Q t X + 1 Q t A 5 - Q T nAx - |Q T nA 2 5 


QX + ^QA5 - QnAx - jQnA 2 8 
I 3 (l - 2nA + n 2 A 2 )5 + T 1 O 3 


(13. a) 


(13 .b) 


Examination of Eq. (12.k) reveals that x = X' A ; therefore, X a= X' a2 - Neglecting all expressions 


containing A 2 in Eq. (13.b) yields 


P 2 H 


i 4 (y+Pq) Q(x+i A §) 
Q T (X + i A 8) I 3 [8(l-2nA) + qJ 


(13. c) 
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or 


p 2 (-) 



Qb 

I 3 c 


(13.d) 


Since P 2 (-) and P,(-) (see Eq. 12. C) have the same form, and since H and R are constant 


matrices, K 2 has the same form as K, and P 2 (+) has the same form as P,(+) . This rule prevails 
in all future gain and covariance matrices, consequently 


K = 


a n^4 


QlP»Aj 


(14. a) 


The approximations made in deriving the expression for K n have been proven to be justified in 
the numerous runs that were made. As an example, in one run (the data of which will be 
presented in the ensuing) we obtained the following Kalman gain matrix at t 2300 = 23 sec . 


K 


2300 — 


0.9198 

-0.0003 

-0.0004 

- 0.0010 

7.5994 

- 3.1834 

2.3294 


-0.0003 

0.9199 

-0.0004 

-0.0009 

3.1834 

7.5994 

- 2.4496 


-0.0004 

-0.0004 

0.9196 

- 0.0012 

- 2.3294 

2.4496 

7.5994 


- 0.0010 

-0.0009 

- 0.0012 

0.9172 

- 2.4496 

- 2.3294 

- 3.1834 


^,2300 

Kw,2300 


(14.b) 


In this work we consider the continuous time case. In practice, though, measurements are done in 
discrete time. The question is whether the assumption that the time interval, A , is small, can be 
justified in the case where the time between measurements is large. The matrices P k (+) and 

P k+I (-) have the form presented in Eq. (13.d) where a, b, and c are scalars. Therefore, if one 

divides the time interval between two discrete measurements into sub-intervals, one can compute 
the propagation of the covariance matrix between the measurement time points as successive 
covariance propagation through the sub-intervals where the format of Eq. (13.d) is preserved. In 


- 10 - 


doing so, one can always choose A , the sub-interval size, to be as small as desired. Thus it can be 
always chosen such that assumption A 2 « 0 is certainly true. Hence the assumption, on which the 
form of Eq. (14. a) was obtained, is fully justified. In conclusion, even propagation over a large 
time interval results in a covariance matrix which has the format of Eq. (13.d). Therefore the 
Kalman gain matrix takes the form of Eq. (14.b) even when the time interval between two 
consecutive measurements is large. 

The realization that K q and K 0) take the special form discussed hitherto, has a far 

reaching consequence. Because of this special form, one needs to tune at most five scalars when 
running the filter of Eq. (9). One scalar has to be determined for each one of the gains. The other 
three are needed if the Markov process in Eq. (9) may have three different components. In such 
case the matrix nl 3 takes the form diag(n, n 2 n 3 ) and then the three n i rather than a single n 

have to be evaluated. In many cases though, K q is simply I 4 and n has indeed a single value, 

thus only two parameters have to be tuned, which makes the filter very attractive. Moreover, 
since the gain is computed directly, there is no need to compute the recursive covariance part of 
the Kalman Filter algorithm if the covariance is evaluated to solely compute the Kalman gain. 

The general results K q = al 4 and K u = Q^p n A can easily be reasoned as follows. Consider 
the upper part of Eq. (9); namely 

4 = jQd> + K q (q m -q) (15) 

When at some point 6) and q equal, respectively, to and q m , then the solution of the equation 
q = 1/2 Qd) settles on q m (and it is assumed that the latter is very accurate). However, when this 
is not the case, the computation of q is driven towards the correct value by the difference 
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(q m - q) . Since this difference directly affects the computation of the estimated quaternion, it is 
logical that each one of the 4 components, and no other component of the difference, drive the 
corresponding component of q ; that is, K q is a diagonal matrix. Moreover, as there is no reason 
to prefer any component over the other 3, K q has to be of the form K q = al 4 where a is some 
constant. In particular, a can be equal to 1 ; that is, K q = I 4 . 

Let us consider now the lower part of Eq. (9); namely 

= -Na> + K {0 (q m -q) (16) 

As will be shown later (Eq. 1 8.c), a least square estimate of d> is given by 

cb = 2Q T q (17) 

When to is equal to to and q = q m , q will stay equal to q m (see Eq. 15); however, as indicated 
in Eq. (16), when q * q m , the difference between the two serves as a signal for changing to in 
the correct direction. According to Eq. (17) the rate of change of a quaternion influences the 
estimate of to through 2Q T . This corresponds to the conclusion expressed in Eq. (14.a) that K 0) 
is proportional to Q T . 

In the next section, where we compute the angular rate using a simple feedback loop, we 
will take advantage of the approximations of K q and K o) . 

PSELIKA BASED RATE EXTRACTION FEEDBACK LOOP 

Following the discussion on the form of the gains K q and K m in the preceding section, 

we use in the block diagram of Fig. 1 the appropriate gain forms. The resultant block diagram is 
shown in Fig. 2. Since q is very accurate, we use Q(q m )that we denote by Q ra rather than 
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Q(q) . Also note that nl 3 = N . Only two parameters; namely, u and nneed to be determined by 
tuning. 



Fig. 2: Block Diagram Representation of the PSELIKA Filter Rate Estimator 


We ran a simulation with o = 5500 and n = 0.1 1/sec. We chose a baseline angular 
velocity pattern that was composed of various features. It included an oscillatory part, abrupt 
changes, ramps and a straight line, all at different levels. This angular velocity generated perfect 

q m 


True and Estimated Rate 



Time (sec) 

Fig. 3: True (red line) and Estimated Rate History for the PSELIKA Derived Loop 
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Fig. 3 (red line) 11 presents the nominal baseline angular velocity components together with their 
estimates. The components of the rate estimation error are shown in Fig. 4. The difference 
between the measured and estimated quaternions is shown in Fig. 5. 


Rate Estimation Error 
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Fig. 4: Rate Estimation Error for the PSELIKA Derived Loop 


SIMPLE PASSIVE FEEDBACK SYSTEMS FOR ANGULAR RATE EXTRACTION 

Borrowing the notion that the angular rate is estimated by a KF in a feedback manner, one 
can raise the question whether the rate can be estimated using a simple feedback loop when the 
attitude is available almost continuously. The logic behind this concept is as follows. If a 
feedback loop contains a node which is an input to an integrator, then the variable at the node is 

1 This rate will serve as a baseline rate for all the tests that follow. 
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the derivative of the output of the integrator. If that output is subtracted from the measured 
quaternion, and the resulting difference is fed forward into the node in a way that drives the 
difference to zero, or almost zero, then the output of the node is approximately equal to the 
measured quaternion. Therefore the variable at the node is very close to the quaternion time 
derivative. When the latter is known, then a good estimate of the angular rate can be computed. 
Before proceeding with the exposition of the feedback control loop idea, we first show how the 
angular rate can be computed from the quaternion time derivative. 



Fig. 5: Quaternion Estimation Error for the PSELIKA Derived Loop 


Angular Rate Determination from Quaternion Rate 

Eq. (2.a) is repeated here 

q=iQ« (2. a) 

where Q is as given in Eq. (2.b). If q and q are known, then we can obtain the following least 

square estimate 14 of to 
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d> = 2Q # q 


(18.a) 


where Q # is the pseudo-inverse of Q , given by 

Q # =(Q t Q)-'Q t (18. b) 

Using the fact that Q T Q = I 3 and Eqs. (18), one obtains the following least square estimate of the 
rate vector 

d> = 2Q T q (18. c) 

We note that this form corresponds to form of the K m gain of Eqs. (14). 

The use of Eq. (18.c) requires the differentiation of q. However, as mentioned in the 
Introduction section, this introduces unwarranted noise in the computed angular rate. However, 
when used in conjunction with the quaternion derived through feedback, the differentiation is 
eliminated. The elementary feedback method for estimating the angular rate is shown next. 

A Simple Gain Feedback Loop 

The ability to achieve satisfactory rate estimates using a relatively simple feedback loop 
with a time varying gain such as in Fig. 2, raises the question whether it is possible to perform the 
task with an even simpler feedback loop. The simplest control loop for deriving q from the 

measured quaternion, q m , and computing <h from it, is depicted in Fig. 6. As mentioned 
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earlier, when q follows q m , q follows q , thus the estimation of to is possible, as indicated in 

Eq. (18.c). When K is a diagonal matrix, the loop is stable for any positive values on the main 
diagonal of K. A simulation was run of the feedback loop in Fig. 6. The gain matrix, K , in this 
simulation was 7 • I 4 . The components of the true and estimated rates are presented in Fig. 7 and 

the difference between them is shown in Fig. 8. Figure 9 presents the quaternion estimation error 
components. 
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Fig. 7: True (red line) and Estimated Rate History When Using the Simple Feedback loop 


An Integral Feedback Loop 

The advantage of the previous scheme is its simplicity. Its disadvantage is its limited 
accuracy. This stems from the fact that for very high gains s is a very small signal that contains 
digital and sensor noise. This noise is then multiplied by a large gain to produce a very noisy q . 
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As a result, the estimated rate is of low fidelity and is very noisy too. On the other hand, for very 
low gains, s is not small; that is, q does not follow q m closely, hence q is not close enough to 
q m ; therefore, the estimate of the angular velocity is of low quality. If, however, we can isolate 
q from £ , then we can keep e very low, thereby force q to follow q m very closely, and yet 

maintain the closeness of q to q m . This can be done by replacing the pure gain in the loop of 

Fig. 6 by a gain and an integrator; that is, replacing K by K/s. This, however, will generate a 
marginally stable loop with poles on the imaginary axis. To stabilize this loop, we can shift the 
marginally stable poles to the left by changing K/s to K/(s + a) which is, basically, a low-pass 
filter. This loop block diagram is shown in Fig. 10. One can tune Kand a to achieve the 
desirable performance. The components of the rate and quaternion estimation errors for 
K = 10,000 and a = 100 are shown, respectively, in Figs. 1 1 and 12. A comparison between the 
results of Fig. 12 and the results presented in Fig. 8 indicates that the added low-pass filter 
reduces the rate estimation error and eliminates the numerical noise as well. 



Fig. 10: A Stable 2 nd Order Feedback Loop for Computing a Rate 

Estimate. 
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A comparison between Figs. 9 and 13 reveals that £ is about two orders of magnitude smaller 
when applying the feedback integral-loop. Thus indeed, we succeeded in forcing q to follow 

q m very closely, while maintaining the closeness of q to q m . As a consequence the rate 
estimation error is much smaller, as indicated by the comparison of Figs. 8 and 12. 


True and Estimated Rate 



0 50 100 150 200 


Time (sec) 

Fig. 11: True (red line) and Estimated Rate History When Using the Integrated Feedback 


RATE-ESTIMATION WITH NOISY QUATERNION MEASUREMENTS 

In the previous sections we considered ideal quaternion measurements. It is interesting to 
examine realistic cases where the measurements include measurement errors and learn of their 
effects on the rate estimation techniques presented hitherto. A typical lateral measurement error 
of an AST is about 40 arcsec. This translates to about 0.0002 rad (see the Appendix). Such white- 
noise error at a rate of 1 sec translates to an equivalent white noise error of 0.00001 rad when the 
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dwz (deg/sec) dwy (deg/sec) dwx (deg/sec) 









measurements are taken every 0.01 sec, which is the step size of our simulation runs (see Eq. 8.3- 
24 in Ref. 14). Using this value as the standard deviation of numbers drawn from a Gaussian 
distributed random number generation; we reran the cases presented in the previous sections. In 
Figs. 14 and 15 we present the true and the estimated rates when we ran the PSELIKA rate 
estimation loop shown in Fig. 2. Similarly, Figs. 16 and 17 present the same for the simple gain 
feedback loop of fig. 6, and finally, the plots in Figs. 18 and 19 present the results for the stable 
second order feedback loop shown in Fig. 10. 

Adding noise to the measurements induces noise in the estimated rates, but the estimated 
rate components closely follow the true rate components. 


True and Estimated Rate 
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Fig. 14: True (red line) and Estimated Rate History for the PSELIKA Derived Loop 
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Rate Estimation Error 



Fig. 19: Rate Estimation Error When Using the Integrated Feedback 

CONCLUSIONS 

In this work we investigated the possibility of extracting the angular velocity vector from 
quaternion measurements without resorting to Kalman filtering. Avoiding Kalman filtering rids 
us of the cumbersome recursive computation of the covariance matrix. It was found 
experimentally, and then justified, that when estimating the angular velocity from quaternion 
measurements using the Pseudo-Linear Kalman filter and a simple dynamics model, the Kalman 
gain has a peculiar structure. The part of the gain which influences the quaternion estimate is 
proportional to, and usually very close to, the four dimensional identity matrix. The part that 
influences the rate estimate is proportional to the transpose of the Q matrix from the kinematics 
quaternion equation. Using this quality, a passive feedback loop was designed in which only two 
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scalar parameters had to be tuned in order to obtain good rate estimates. As mentioned, this 
feedback loop uses the Q matrix. Therefore, although it is a passive system, it has a time-varying 
gain, and thus it is a linear system which is time varying. 

Two other feedback loops were also introduced which are passive, as well as constant 
parameter systems. The systems consist of four simple decoupled loops. The first feedback 
system employs a simple diagonal gain matrix which multiplies the four components of the 
difference between the estimated and measured quaternion. Each element of the amplified 
difference is fed into an integrator, the output of which is a component of the estimated 
quaternion. Therefore, the input to each integrator is the derivative of the estimated quaternion. 
With these derivatives on hand, the rate estimate is obtained by a known simple linear, time- 
varying transformation. The second feedback system constitutes an improvement on the first one. 
It contains an added pole that enables a more accurate estimation of the angular rate. 

Simulation results indicate that all three feedback systems presented in this work are adequate 
for deriving the angular rate vector from frequent quaternion measurements for the purpose of 
attitude control loop damping. 

Although we only considered quaternion measurements as input to the rate estimation 
feedback loops, vector measurements can be handled too. This is done by applying the method 
presented in reference 15 where a measured unit vector and its expression in the reference 
coordinate system are converted to a pseudo-measured quaternion. 
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Appendix 


The purpose of this appendix is to show how we determined the standard deviation of the 
quaternion measurement error when given the AST angular error. Let the pointing error of its 
Euler vector be d(p x , dcp y , dcp z (note that these errors are not independent), and let the angular 

error about this vector be denoted by d(p. The equations relating the nominal quaternion 
components to the nominal components of the Euler vector are 


• 9 

q, = sm-^--cosa x 

• 9 

q 2 = sin-j-cosa y 


q 3 = sin— -cos a z 


q 4 = cos 


9 

2 

9 


(A.l) 


where cos(a 1 ) = dcp s / cp i = x, y, z . Perturbation of the quaternion components yields 


dq, =^cos~cosa x -dcp-siny -sina x da x 

dq, = — cos— -cosa v •dcp-sin— -sina v -da v 
H2 2 2 y 2 y y 

dq, =— cos— • cosol - dcp-sin— -sin a 7 - da, 

M3 2 2 2 

dq 4 = -—sin— -dtp 
44 2 2 


(A.2) 


Assuming worst case for each component (sine as well as cosine functions equal 1 and the 
negative term eliminated) we obtain: 
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(A.3) 


dq, = — dcp 
d q 2 =^(p 
dq 3 =|dcp 

dq 4 =_ ~dcp 

We were given the error value d 9 = 40 arcsec which translates to about 2 e-4 rad. We 
conservatively choose the standard deviation of each one of the components to be equal to the 
differential of that component, hence 

CT i,n = ~d(p = ~2 e-4 rad = 1 e-4 rad i = 1,2, 3, 4 (A.4) 

where n is the discrete time step, and, as before, i denotes the quaternion component number. 
Then 


Q jn =a 2 = le-8rad 2 i=l,2, 3, 4 (A.6) 

Using the approximate formula for transforming variance (discrete sequence) to power spectral 
density (continuous process) (Ref. 14, Eq. 8.3-24) we obtain 

Qi.n = Qi • At n (A.7) 

For a time interval between two consecutive measurements of 1 sec we obtain 


q. = =q =1 e -8rad 2 (A.8) 

At n 1 

Using the same transformation formula as before, we find that for a measurement time interval of 
0.01 sec, the equivalent variance is 

Qi „ = Qi • 0.01 = 1 e-10 rad 2 (A.9) 
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hence 


Oj,„ = Qi/ 2 = 1 e-5 rad (A. 10) 

This is the standard deviation that we use to draw numbers for simulating the quaternion noise. 
We draw the numbers from a Gaussian random number generator of zero mean. 
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